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ABSTRACT 

The stationary hydrodynamic equations for the transonic accretion disks and flows 
around rotating black holes are presented by using the Kerr-Schild coordinate where 
there is no coordinate singularity at the event horizon. We use two types of the causal 
viscosity prescription, and the boundary conditions for the transonic accretion flows 
are given at the sonic point. For one type of the causal viscosity prescription we also 
add the boundary conditions at the viscous point where the accreting radial velocity is 
nearly equal to the viscous diffusion velocity. Based on the formalism for the transonic 
accretion disks, after we present the calculation method of the transonic solutions, 
the horizon-penetrating transonic solutions which smoothly pass the event horizon 
are calculated for several types of the accretion flow models: the ideal isothermal 
flows, the ideal and the viscous polytropic flows, the advection dominated accretion 
flows (ADAFs) with the relativistic equation of state, the adiabatic accretion disks, the 
standard accretion disks, the supercritical accretion disks. These solutions are obtained 
for both non-rotating and rotating black holes. The calculated accretion flows plunge 
into black hole with finite three velocity smaller than the speed of light even at the 
event horizon or inside the horizon, and the angular velocities of the accretion flow at 
the horizon are generally different from the angular velocity of the frame-dragging due 
to the black hole's rotation. These features contrast to the results obtained by using 
the Boyer-Lindquist coordinate with the coordinate singularity at the horizon. 
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1 INTRODUCTION 

Elucidating the nature of the strong-gravity region around a black hole is one of the greatest challenges in astrophysics in 
this century. The most of the gravitational energy is converted into another forms of energies such as the kinetic, the thermal 
and/or the radiation energy, and especially the emission from the vicinity of the event horizon contains information about 
physics of the strong-gravity region such as the physical parameters of a black hole. Since such emission by e.g. photons or 
neutrinos is usually produced in the accretion flows plunging into the black hole, the structure of the accretion flows in the 
vicinity of the black hole's horizon is frequently required to be solved precisely. Especially, the effects of the frame dragging 
due to the black hole's rotation are taken into account only when the fully general relativistic calculations are done. 

In the past studies, for many astrophysical systems and situations, the stati onary solutions for the structure of the 



transonic accretion flows near the horizon are solved for the standard accretion disks dNovikoy fc Thorn j| 19731: IPaee fc T: 



1974). the advection-dominated accretion flows (IChakrabarti|[l996l: lAbramowicz et aLlll997l : Ijaroszvski fc Kurpiewskil | l997l : 
Gammie fc Popham 19981 : Poph am fc Gammie"l998';'Manmoto|l2000l)_^ t he po lytropic accretion flows (jPeitz fc App]||l997h ~the 
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super-critical accretion flows Ifioloborodoy 1998 ; Shimura fc Manmotol|2003l '). and the hypercritical accretion flow model for 



the neutrino-dominated accretion flow (jPopham. Wooslev fc Frvei 19991 ). All these works use the Boyer-Lindquist coordinate 



where the coordinate singularity exists at the horizon. Due to the coordinate singularity, some physical values of the accretion 
flow based on the calculations using the Boyer-Lindquist coordinate are not realistic. For example, the radial component of 
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three velocity of the accretion flow equals to the speed of light at the horizon, and the corresponding gamma factor diverges 
just outside the horizon. But these features are not the cases for the realistic accretion flow, i.e., the accretion flows plunge 
into the horizon with some finite velocity smaller than the speed of light, because from the point of view of the local observer 
moving along the fluid motion, the horizon is not a special location. So, one of the possible natural next step is to extend 
these past studies to the formulation using the Kerr-Schild coordinate where there is no coordinate singularity at the horizon. 
Our calculations in this study by using the Kerr-Schild coordinate avoiding the coordinate singularity show the radial velocity 
smaller than the speed of light at the horizon. In this study, we first show all the explicit formulation for the transonic accretion 
fiows written by the Kerr-Schild coordinate, and then calculate the horizon-penetrating transonic solutions for the accretion 
fiow which do not have any singularity at the horizon. 

The Kerr-Schild coordinates are frequently used in the past study, e specially for the dynamical numerical calculations of 

the hydrodynamics or the magnetohydrodynamics around the black holes (jPapadopoulos fc Fontll99^: Font. Ibanez fc Papadopoulosh 



tne nyaroaynamics or tne magnetonyaroayna,mics aroung tne piacK notes (|-rapaQOPOuios ez -fontiiij^aiMront. ipanez &: rapaaopouiosi 
1994lcookl2000l : lFontl200(]| : iKomissaro^^boOll : ICammie. McKinnev fc T6thll2003l : lKomissarovll2004l : lGammie. Shapiro fc McKinnev^ 



2004). While these studies m ainly concentrate on the dynamical calculations of the accretion flows around the black hole 



Papadopoulos fc Font j 19981 ) also give the simple solutions for the stationary accretion fiows plunging into the black hole 



which exhibit the maximum value of the radial speed at the horizon. In our calculations, we found that this feature is not a 
general statement in the viscous accretion fiows, i.e., the maximum value of the radial velocity is not generally achieved at 
the position of the horizon. 

In the present study, we concentrate on the stationary axisymmetric accretion fiow in the equatorial plane. In this 
study, we use two types of the causal viscosity prescription. One is the simple treatment of the kinematic viscosity so as to 
vanish the shear stress at the horizon. The other is based on the description of the shear stress measured in the fiuid's rest 



frame. The latter ty p e of t he causal viscosity is considered bv lPapaloizou fc Szuszkiewicd (|l994l l for the Newtonian case and 



Gammie fc Popham ( 19981 ) for the relativistic case. By using the latter types of the causal viscosi ty prescription, we do not 



need to put the boundary condition on the horizon, such as zero-boundary condition used in, e.g.. iNaravan. Kato fc Honma 



(1997]). In this study, we give the formulation with the effects of the heat fiux (or the radiation cooling), the heat inertia and 
the relativistic enthalpy, and the relatively general forms of the equation of state. Based on this formalism, we give the sample 
transonic solutions for three types of accretion fiow models: ideal isothermal disks, polytropic disks and advection-dominated 
accretion fiows (ADAFs) with the relativistic equation of state, the adiabatic accretion disks, the standard accretion disks, 
and the supercritical accretion disks. 

We give the preliminaries for the calculations of the transonic solutions by using the Kerr-Schild coordinate in §2 containing 
such as the background metric and the frame and the frame transformation used in this study. In |3l the basic equations 
are given; the mass conservation f ^3.ip . the radial momentum equation fi i3.2p . the angular momentum equation (i 33.3p . the 
equation for vertical structure which is the momentum conservation in 6'-direction ('i i3.4p . the thermodynamic equation and 
the energy equation (i ]3.5|l . the turbulent shear stress based on the causal viscosity (i ]3.6|l . The boundary conditions for the 
viscous transonic solutions are given in 52] the boundary conditions at the sonic point (JSTT)) and the viscous point ( ij4.2p . 
In !j5] we summarize the coupled differential equations to be solved. We give the calculation procedures in ^ Formulation 
and/or Numerical solutions for the horizon-penetrating solutions of the transonic accretion flows for different types of the 
accretion flow models are presented in SJT] the formulation and the numerical solutions for the ideal isothermal accretion 
flow ( ii7.ip . the polytropic disks (i )7.2p . the ADAF with relativistic equation of state fi l7.3p . the adiabatic accretion disk and 
the standard accretion disk fi i7.4p . the formulation for the supercritical accretion disk ( ^7.5p . We give concluding remarks 
in the last section. While the basic structure of the basic equations for the transonic accretion flows are simple, the explicit 
expressions for some formula are lengthy. In order to clearly see the outline of the calculations, we put the details of lengthy 
formula in the Appendix and in the main body only the important formula are given. 

2 METRIC, REFERENCE FRAME AND VELOCITY FIELDS 
2.1 Background Metric 

Throughout the present study, we assume the background geometry around the rotating black hole written by the Kerr-Schild 
coordinate described as 

ds^ = -a'^dt^ + -t^3{dx' + l3'dt){dx^ + /3^dt), (1) 

where i, j — r, 8, (j), and the nonzero components of the lapse function a, the shift vector /3* and the spatial matrix 7ij are 
given in the geometric units as 



_ / 2mr\-i''2 2mr/S 

' ' ^ l + 2mr/E' 



2mr ^ A sin"' 61 2„/ 2mr\ , . 

7rr = 1 4- ^-,799 = E, 700 = ^ , 7^0 = 70r = -« sm g (1 -|- 1 . (2) 
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Here, we use the geometric mass m = GM/c^, E = + cos^ 6, A ^ — 2Mr + a? and A = (r^ + a^)'^ — a^Asin^ 0, where 
M is the black hole mass, G is the gravitational constant and c is the speed of light. Explicit forms of nonzero components 
of metric g^v and its inverse g'^^ are calculated in Appendix|X] The position of the outer and inner horizon, r±, is calculated 
from A = as r± = m±{m? —a?Y^^. The angular velocity of the frame dragging due to the black hole's rotation is calculated 
as = ^gt(t>/ g4>(t> ~ 2mar/A. Although some past studies use the metric in the equatorial plane, i.e. 6 — it/2, we basically 
formulate the basic equations by using the metric including 9. This is because in the calculations of the vertical structure of 
the accretion disk which are performed in Sec. 13.41 and Appendix [F] need the metric including 9, and we would like to have 
the consistency of the notation of the metric throughout the paper. But, we actually calculate the transonic solutions of the 
accretion flow in Sec[7]by setting 9 = n/2. 



2.2 Reference Frames 

In this study, we use three types of reference frames. The first is the Kerr-Schild coordinate frame (KSF), in which most of 
our calculations are done. The second is the fluid's rest frame (FRF), an orthonormal tetrad basis carried by observers moving 
along the fluid. The components of the four velocity measured in the FRF are described as 



u 



= 1, u^"' = U(i) =0, (z = r, 9, cj)) (3) 



where the bracket denote the physical quantities measured in the FRF. The third frames are calculated from a stationary 
congruences formed by observers with a future-directed unit vector orthogonal to t =constant whose components are given as 

ut = -a, Ui =0, (i = r, 9, 0), (4) 

and 

y.-^, u' = -a~^l3', (i^r,9,(P), (5) 



It 



respectively. For this congruences, since the vorticity tensor vanishes [e.g., A. 10. 2 in lFrolov fc Novikovl (|1998| )]. this congruences 
is the congruences of locally non-rotating observers. So, this frame is usually called as the locally non-rotating reference frame 
(LNRF). By using the Boyer-Lindquist coordinate, this observer is moving with the angular velocity of the frame-dragging due 
to the black hole's rotation (Bardccn 1970; Bardecn, Press & Toukolsky 1972). On the other hand, by using the Kerr-Schild 
coordinate, since /J*" 7^ and f3^ — fJ''' — 0, the observers with = — is radially falling with — u''' — 0. For such 
observers, the nonzero components of the covariant and the contravariant four velocities are given as 

LNRF t 1 r fa\ 

Q a 

We can easily show that the LNRF is an orthonormal tetrad basis carried by the observer moving with u^j^j^p. The physical 
quantities measured in the LNRF are described by using the hat such as w'', u^, etc. In the Kerr-Schild coordinate, since the 
congruences of the observers moving with the angular velocity of the frame-dragging have the singularity at the event horizon 
as shown in the Appendix [B] we do not use such congruences in our calculations. 



2.3 Frame Transformations and Velocity Fields 

The physical quantities measured in the KSF are transformed to those in the LNRF by tetrads e^^ and e**-. For example, the 
four velocity is transformed as u'^ — e^-u" and — ei^ut,. The explicit expressions of the tetrad components e^j" and e''- are 
given in Appendix [Cl 

The FRF usually moves with some radial and azimuthal velocities with respect to the LNRF. We newly deflned the radial 
velocity Vr and the rotational velocity such that the FRF moves with the radial velocity Vr and the azimuthal velocity 
ii^ with respect to the LNRF. By using these velocities, the physical quantities in the LNRF are transformed to those in the 
FRF by two-dimensional Lorentz transformation e-'^' and e"(\) with the radial velocity Vr and the azimuthal velocity i)^. 
Here, e-'^' and e"(\) are the transformation matrices denoting the two-dimensional Lorentz transformations, and the explicit 
expressions of these matrices are also given in Appendix [Cl 

By using the tetrads described in Appendix [Cj all the covariant and contravariant components of the four velocity in 
the KSF are calculated as — e''^e" j^jii'"' and = e/j"e^''''''M(^), and described by using the radial velocity Vr and the 
azimuthal velocity of the FRF measured in the LNRF as shown in Appendix |D] Inversely, the radial velocity Vr and the 
azimuthal velocity are described by the four velocities, u'^ and Ufj,, measured in the KSF as 

where 7 = (1 — — i'"^)"^^'^ = olu^ . Since Vr is the radial velocity measured in the LNRF which is radially falling with iilnrfi 
the radial velocity Vr can generally have both the positive and the negative values for the radially falling accretion flows. 
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2.4 u* and 

When the transonic solutions are calculated later, we solve the differential equations for the radial four velocity and 
the angular momentum £{= u^). Therefore, it is convenient to express the angular velocity fl and it* which are frequently 
used in the formula in the following sections by and £. From the normalization of the four velocity, u^^Uf^ — —1, we can 
obtain the quadratic equation of u' as Ao{u^)^ + 2i3ow* + Co = where Ao — —c? + {fy^ 1^^'^ , Bo — {ff /'y^'^)u'^ and 
Co = {vI'Y /^^^ + ^^/7^0 + 1. From the quadratic equation, we have the solution for u* > as 

t C'o 



(8) 



where Do ~ Bq — AoCo- In this study, we consider the accretion flow with < 0. For such flows. Bo < 0, and then u* > 0. 
We can also show that Do > for the region r > r+, and for the region r < rg, Dq > only when £^ < [{u^)'^ + 1/T,]A/{~A). 
From u* calculated above, the angular velocity Q is calculated as 

^^^^e_Jr^ (9) 

which is derived from £ = y^^ii'^. Here, is calculated from Eq. (|8}. We can also calculate Ur and £{= —ut) from and £ 
as Ur — [gtr + gr4,^)u^ + QrvU^ and £ — [gtt + gt4,^)u^ + gtrvJ' where w' and are calculated by Eqs. (|8]) and ((Ojl. 

2.5 Transformation of Four Velocities written by Kerr-Schild Coordinate and Boyer-Lindquist Coordinate 

The transformation of the four velocities calculated by using the Boyer-Lindquist coordinate and the Kerr-Schild coordinate 
are given by 

t t 2mr rrreethsO'r 
^^BL = "kS a""'^^' '^BL^ftKS, ^BL^ftKS, ^BL = ""kS ~ ^^KS , (10) 

BL KS BL KS , 2mr KS , " KS BL KS BL KS /-.-.n 

Ut =Ut , Ur = Ur + —^^t + -^^4, , "U-S ^ U0 , U^ ^ U^ . (11) 

Here, "BL" and "KS" denote the physical quantities calculated by using the Boyer-Lindquist coordinate and the Kerr-Schild 
coordinate, respectively. 



3 BASIC EQUATIONS 

The basic equations for the relativistic hydrodynamics are the baryon-mass conservation (pou'');^ — and the energy- 
momentum conservation T.^t" = 0, where po is the rest-mass density and T'"' is the energy-momentum tensor. Dynamical 
basic equations except the baryon mass conservation are calculated from the energy-momentum tensor, y". We use the 
energy-momentum tensor written as, 

= poriuu + pg +v +q u +q w, (12) 

where p is the pressure, 77 = (po + ^'+p)/po is the relativistic enthalpy, u is the internal energy, t^" is the viscous stress-energy 
tensor and is the heat-flux four vector. In the present study, we do not include the heat flux term in the energy-momentum 
tensor. 

One of the natural form of the shear stres s, f^*^ , is the Navier-Stokes shear stress. The relativistic Navier-Stokes shear 



stress, f", is written as ( Misner. Thorne fc Wheeler 1973f) 



t"" = ~2\a'"' - CQh"'' , (13) 

^ The transformation law given by Eq. mOI I is calculated from the lapse function a, the shift vector /3'' (or and the matrix 7^^- (or 
7'^ ) for the Kerr-Schild coordinate as 

t _ t r <!> _ ■t> , f ^r,t,\ a'^Y'' r 

"bL-"kS Q,2^rr _ (^r)2"KS' «BL ~ "^S + I — ] „2^rr „ lar\2'^KS' 



which are derived from the metric expressed as 



^rr 



dt : — TTrdr 

a2^rr _ (^r)2 



a2^rr _ (•^r)2 
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where A is the coefficient of dynamic viscosity, ^ is the coefficient of bulk viscosity, /i^" = g*"^ + u'^u'^ is the projection tensor. 



e ; 



is the expansion of the ffuid world line, and a'"' is the shear rate of the ffuid which is calculated as 



7 a I 



(14) 
(15) 



where a„ 



is the four acceleration. In this study, w e do not take the shear stres s wr itten by this form. In s tead, we use 
the Kerr-Schild coordinate version of the shear stress used in ICammie fc PophamI and lPopham fc Gammi^ l|l99^ which 
allows angular momentum transport and preserve causality. We evaluate the shear stress in the FRF and assume that all the 
components of the shear stress except t{r){<p) ~ ^{<t>){r) are null in the FRF. Based on this assumption, the shear stress measured 



1 t 



{r)(4>)- 



in the KSF is calculated by using the tetrads connecting the KSF and the FRF, e.g.. 
The explicit forms of the shear stress in the FRF used in this study is given in Sec. 

In this study, we consider the stationary, axisymmetric and equatorially symmetric global accretion ffow in the equatorial 
plane, i.e., we assume — 0. We also assume that the effects of the bulk viscosity is negligible. In the following sections, 
we derive the basic equatio ns written by the Kerr-Schild coordinate by using the vertical averaging procedures used in, e.g.. 



Gammie fc PophamI (|l998l '). around the equatorial plane. 



3.1 Mass Conservation and Mass-Energy Flux 

The equation for the baryon mass conservation is written as 

(Pou'');M = 0, 



(16) 



where po is the rest-mass density and u'^ is the four velocity. By averaging the physical quantities around the equatorial plane, 
the mass-accretion rate M is calculated as 



(17) 



where Hg is the half-thickness of the accretion flow in the ^-direction which is calculated in Sec. 13.41 and ^/—g = r^. When we 
calculate the global structure of the accretion flow, we normalize the rest-mass density, po, by setting A/ — 1, i.e., the mass 
conservation is written as 



Anr^ Hgpou'^ = 1. 



(18) 



From the projection of the energy-momentum conservation, r'?J^ = 0, onto t-component, i.e., hl^T'^^'^ — 0, with the vertical 
averaging calculations, we get 



^nHg. 



M 



-(t t+u qt-q £) 



^TvHgy 



M 



-£q 



From Eq. (|19|l . we obtain 

A-KHg^ 



r]£ ~ eo + 



M 



'-t' 



Qe, 



where eo is the specific energy of the fiow and Qe represents the effects of the heat flux defined as 



Qe 



ATvHg^/^( 



■£q dr ^ 



M M 
In the case of no heat flux, Eq. (|19|l is reduced to 
E = Mr]£ + A-Kr'^Hgft, 



{uqt - q''£) . 



(19) 



(20) 



(21) 



(22) 



where E is the mass-energy flux corresponding to the rate of change of the black hole mass if measured at the horizon. With 
the mass conservation, we obtain the specific energy eo as 



E_ 

M 



r]£ + 



t t 
poW 



(23) 



where eo is the specific energy of the accreting matter. When the velocity of the accreting matter is non-relativistic and cold 
i.e. rj — 1 where the thermal energy of the matter is much lower than the rest-mass energy, the specific energy become unity, 
i.e. E ~ M. 
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3.2 Radial Momentum Conservation 

The equation for the radial momentum conservation is obtained by the projection of the equation for the energy momentum 
conservation, T'^J^ = 0, into r-direction, i.e., h^T^]^ — 0. We can write down h^^T^^ — as, 

a = -f+nm, 24 

pot] dr 

where is the radial compo nent of the four-acceleration of the fluid, a'^ = 11^.^,11", and tihi includes the effects of the heat 
inertia which are discussed bv lBeloborodov. Abramowicz fc Novikovl l|l997r ). 

The radial component of the four-acceleration, {= u^.^^u'^), is calculated as 

rdvJ' 

a' = u' — riacc, (25) 

dr 

where we decompose Wacc into three parts as 

Wacc = nKp(r2 — ^'k)(^ — ^k) + '^BL + WKS . (26) 

Here, Q, = /u^ is the angular velocity and Q,^ are the Keplerian angular momentum described as Q,^ — ±m^^'^ / {r^^^±am^^'^) 
which are the solutions of g^^^r^'^ +'2gt(i>,r^ + Qtt.r = 0. The term including ukp measures the deviation of the angular velocity 
from the Keplerian angular velocity. The terms wkp, tibl and uks are given as 

riKP = ^g''''94>4;riuf, (27) 

nBh = -is"ffrr,r(M'')^, (28) 
TlKS = — (ff"'5tr,r + S"^gr</.,r)(lt'^)^ — [g*"" (ptt ,r + fft ,^>,rf^) + g'''* (fft ,^>,r + 50</),r fi)] u'u'' . 

(29) 

Since uks contains <;*' and g^"^ which are null for the Kerr metric written by the Boyer-Lindquist coordinate, this term is 
newly calculated term in this study which use the Kerr-Schild coordinate. On the other hand, the general form of nm is 
described as 



riHi = 

PoV 
1 

poV 



-[Kt^:: + h''^{q^U- +q-U^).,.] (30) 



<P — -(St ' + u^u ]u + t -I- q-yU ' + g 6 -|- u.^q 



(31) 



where "I> = — cr^^t'"' is the dissipation function which is calculated in Sec. I3.6l based on the shear stress measured in the FRF, 
and — (l/3)0tXy represents the compressive heating rate. 

If the effects of the dissipation function, <1>, is dominated, nni is reduced to 

r 

11 

nm = $, (32) 

PoV 

and we use this expression in this study. 

Finally, from the equation for the radial momentum conservation, we can derive the equation for du^ /dr as, 

rdu'^ K''' dp 

U — ;— = ; hriacc-l-nHI. (33) 

dr poT] dr 

This equation is used to derive equations which determine the boundary conditions for the sonic point and the viscous point 
in Sec. |4] In this study, we solve the radial component of the four velocity in the KSR instead of the radial velocity Vr 
measured in the LNRF when we solve the transonic solutions. This is because while Vr have the negative and positive values 
as denoted above, is always negative for the accretion flow. Thus, we choose as one of the basic dynamic variables when 
solving the transonic flows. 



3.3 Angular Momentum Conservation 

The equation for the angular momentum conservation is obtained by the projection of the equation for the energy momentum 
conservation, T'f^ — 0, into 0-direction, i.e., h'j^T'^J^ = 0. By using the vertical averaging procedure, we can write down 

htT^fi: ^ as. 



rj£ (t 4, + u q^~q e) 



M 



-iq 



(34) 



From Eq. (|34|l . we obtain 
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-qi-j -rf t^^Qe, (35) 

where j is the specific angular momentum and Qi represents the effects of the heat flux defined as 

7^ £q dr H iu qs - q I) . (36) 

When the angular momentum is not transported by the heat flux, 

Mj = Mr)l - Anr'^Hgt''^, (37) 

where the flrst term in the left-side hand is the total flux of the angular momentum of the fluid, the second term represents the 
amount of the dissipation due to the shear stress, and the right-hand side is the total inward flux of the angular momentum. 
In this study, we assume no angular momentum is transported by the heat flux, and in this case the shear stress tensor t^^ is 
written as 

tl = -pou'-{r,£-j). (38) 



3.4 Vertical Structure 



The equation for the vertical structure is calculated from the equation of the momentum conservation in 0-direction by 
assuming the hydrostatic equilibrium. The equation for the momentum conservation in ^-direction is obtained by the pro- 
jection of the equation for the energy momentum conservation, T'^^ = 0, into ^-direction, i.e., /i^T^^ — 0. Although the 
calculation procedures for the charac teristic angular scale of the accretion flow. He, is basically same as those used in 
Abramowicz. Lanza fc Percivall ( 19971 ). we use several different assumptions. From /iJ^T'f^^ — 0, by neglecting the effects 



due to the heat flux, i.e. = 0, we obtain 
1 dp 



por} 86 



( 



U _^ ^9 dug 

dr 89 



ug ( rdr] g8ri 
r] \ or do 



PoV 



{q"ue + qgu" 
PoV 



(39) 



Here, the terms including the differen tiation of rj and f^g-^ are newly consi dered terms which are not taking into account in 



Abramowicz. Lanza fc Percivall (1997) which 



1. It is note that 



r; — 1 use the equation for Hg derived by Abramowicz. Lanza fc Percivall ( 1997h assuming rf = 1 
We expand the pressure in ^-direction until the order of cos'^ 9 as, 

cos 6*^2" 



Popham fc Gammiel ( 19981 ) which do not 



p{r,9) =po(r) 



Hg 



(40) 



where po{r) is the pressure in the equatorial plane. This expansion is different from the expansion in Abramowicz. Lanza fc Percival^ 
1 19971 ) by the factor 1/2 before {cos9/Hgf. Fr om this, we can calculate dp/d9 — {2cos9/Hg)po{r). From this, the angular 
half-thickness of the disk is calculated as 



Hi 



\ponl 



UgU r)^y 



I g;v 
poV 



{q"ug + qgu" 



rj por) por] 

In order to expand the denominator of Eq. (|50|l until the order of cos 9, we also approximate ug as 

ue{r, 9) = uei{r) cos 9, 

where ugi is 6-component of in the equatorial plane. By using these expansion until the order of cos 6', we obtain 



U Ug^i, — 



I rdug gdug\ 

r -89 ) 

a ( r , gg 2 \ 
is9 (—It it9i,r + g ugij 



UgU ri,v 

n 



Ug 



8ri 



)8ri 



= I u — h u — 



8r 
UgiU 



89 
9 In 77 
8r 



In the same way, we can also calculate T'^^u^u" and —t'^g-^/ipov) until the order of cos9 as, 

^g^U^u" = cose (F^^lt^-u")!, 

tV,i./{pov) = cose [t^^/ipov)]!, 
{q^ug + qgu'').^/{pori) = COS 9 [{q" Ug + qgu").^ / {pori)]i, 



(41) 
(42) 

(43) 

(44) 

(45) 
(46) 
(47) 



where the coefflcients of terms of the order of cos9 are defined as (Fg^ii^it")!, [t"g.,^/{pori)]i and [{q'^ug + qgu'^)-i, / {por])]\. 
Direct calculations leads the explicit forms of {Tg^Uf_iU^)\ as 



8 Rohta Takahashi 



where £l is calculated as 
£t = f- a{£'^ - 1). 



(48) 



(49) 



From Eg. (|49p . I, = I when a — Q or £ = 1. The form of Eq. (|49p is same as the results of lAbramowicz. Lanza &: Percival 
1 19971 ) where the Boyer-Lindquist coordinate is used. This is because the transformations for t and £ between the Boyer- 
Lindquist coordinate and the Kerr-Schild coordinate are given as ^bl = ^ks and £bl = ^ks- In Appendix |Fl we also show the 
direct derivations of Eq. (I49II which is essentially same as the calculations using the Boyer-Lindquist coordinate as shown in 



Abramowicz. Lanza fc Percivall (119971 ). but several points are different. Then, the most general form for Hg is calculated as 



m 



\ poll J 



U U0^r + lie 



dr / 













)- 


po?7 


1 




1 



(50) 



where all physical values such as 



etc. are evaluated at the equatorial plane. In the present stud y, we assume ug 



an d the negligible effect s for th e viscosity and the heat flux. These assumptions are basically same as 



Abramowicz et al 



(1997) and Gammie fc Popham ( 1998h . For the calculations of the transonic solutions in the later sections, we use the angular 
half-thickness of the disk described as, 

Cs 



Hg 



(51) 



3.5 Energy Equation 

The equation for the local energy conservation is obtained from u^T''^^ = as 
r fdu u + pdpo\ + _ , , 

where 

g+s = <i>-lQt\, (53) 
<l7^.d = -q^-l^a^. (54) 



Elere the dissipation function $ is given in Sec. 13.61 and the second term in the right hand side of gvis side represents the 
compressive heating rate. On the other hand, in the right hand side of q~^^, is the mass-energy flux transported out (in) 

to (from) the outside region, and — g'^a^ is the special relativistic correction to — g*";^ due to the heat inertia of the flux and 
represents the effects of the redshift of the flux. Since the left-hand-side of Eq. (|52|) include the change of the entropy, s, as 

ds 

PoTu —=q+^~ q~^^ , (55) 

where poTds/dr represents the advected energy of the accretion flow, this equation represents the energy balance of the 
accretion flows, i.e. (advection cooling) = (viscous heating)- (radiative cooling). For the isothermal flows or the polytropic flows 
calculated in the later sections, we do not use the energy equation given by Eq. (|52p when we solve the transonic solutions for 
these flows. On the other hand, for the general equation of state where the pressure and the internal energy usually are the 
functions of both the rest-mass density po and the temperature T, the energy equation given by Eq. ((52} is required in order 
to solve the transonic solution. As an example of such cases, we solve the transonic solutions for the advection dominated 
accretion flows with the general relativistic equation of state. 



3.6 Treatment of Viscosity limited by Causality 

The viscosity due to the turbulent motion of magnetic field, fiuids and particles such as photons and neutrinos is usually 
considered in the accretion disk. The effect of viscosity is transported to the finite length with the finite viscous timescale, r„, 
and this transportation is limited by the causality. When the fiuid's velocity approach the speed of light as near the horizon, 
it is expected that the viscous transportation become less effective. In this study, we phenomenologically take into account 
the causal viscous effects. The valid treatment of the causal viscosity will be required in future studies. 



3. 6. 1 Type A Causal Viscosity : Simple treatment of kinematic viscosity 

Here, we consider the kinematic viscosity by taking into account the causality. We use the kinematic viscosity coefficient, u, 
in order that the kinematic viscosity vanish on and inside the horizon which is expressed as 
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1^0 fc 





for 
for 



r+ < r, 
r ^ r+, 



(56) 



where r+ is the radius of the horizon and is the kinematic viscosity coefficient wh en the effects of the causality is not 
considered, and fc is a cut-off function described as ( Naravan f 992l : Peitz fc Appll f997 ). 



(57) 



[l-{v/c^ff, for |fi|^c„, 
0, for |«| > c„. 

Here, v = [1 — ^^'^)^/^ where 7 = qm'. These treatments of the kinematic viscosity is similar to those of Peitz fc Appll (1997), 
but several points and explicit expression are different. By using the kinematic viscosity coefficient defined above, we calculate 
the shear viscosity tensor as Navier-Stokes viscosity described as t^i, — —2pQr]i/o^v 



3.6.2 Type B Causal Viscosity : Shear stress measured in fluid's rest frame 



In relativity, the physical meanings are not usually expressed directly in arbitrary fr ames. The FRF is the mo st natural 
place to evaluate the physical processes. We calculate the shear st ress in the FRF as Gammie fc Popham ( 19981 ) by using 
the relativistic version for the causal stress prescription proposed bv lPapaloizou fc Szuszkiewicz We assume the shear 



stress t 



(0)(r) 



5") in the FRF. The other components of the sh ear stress in the FRF except t(r)(0) and i(0){r) are 



assumed to be null. This treatment is same as Gammie fc Popham ( 1998h . By using the tetrads connecting the KSF and the 
FRF, the shear stress t^^ is calculated as. 



t\ = 2 



r(r) (di) I 

e ^ 'e^^' + e 



'"(0)„ (>•) 



(r){0) 



= FS, 



(58) 



where F 



+ e 



)and 



(0) 



(/5) 



finite value of t 



and tetrads in F are calculated by using the LNRF as e"'^' ' = e"^) — e"^e^^^j (a = r, 
- <f), (3 = r, (j>). The explicit forms of tetrad components are given in Appendix [C] For the 

is not null. This feature is contrasted to the shear stress 



erTe,.,'"' {a = 0, p = r. 
the shear stress measured in the KSF, t 



calculated by using the Boyer-Lindquist coordinate as Gammie fc Popha mllll998l) . see Eq. (60) in lGammie fc PophamI (119981 ). 
The equation for the shear stress S is described as ([Gammie fc Popha ml ll998l) 



dS 
dr 



S-So 



(59) 



where r„ is the relaxation timescale of the viscous diffusion and is the equilibrium value of the shear stress. The relaxation 
timescale is related to the propagation speed of the viscous effects c„ as c„ — (u/tv)^^^ where u is the kinematic viscosity. 
The coefficient of the dynamic viscosity A is described by the kinematic viscosity as A = porfv. In this study, the relaxation 
timescale Tv is assumed as Tv = l/il. From these relations, the propagation speed of the viscous effects , Cv, are described by 



the sound speed, Cs, as c„ = a^^^Cs. These treatments are basically same as Gammie fc Pophaml ( 19981 ) 



From the angular momentum equation (|38|) . the shear stress S is calculated as S = — pou^ {rjl— j —Q 1) / F . By differentiating 
this equation by r, we obtain dS/dr and substitute dS/dr to Eq. (|59|) . Then, the shear stress S is calculated as. 



S = 



dr 



+ 



^dlnry 
dr 



1 dQe 
r/ dr 



dlnF dlnHg 

+ 



dr dr 

where dQe/dr is determined by the heat ffux and defined as 



dQi 
dr 



'9 nO 



M 



£q" + 



M 



From Eq. (|34|, dQe/dr is also calculated as 



dQi 
dr 



r)l- 



M 



-t' 



The equilibrium value of the shear stress 5o is assumed to be the Navier-Stokes value as 
So = -2po77i/cr(,.)(<^). 



(60) 



(61) 



(62) 



(63) 



The shear rate (^(r)(cf>) in the FRF is calculated by using the shear tensor cr^i/ in KSF as cr{r)(</>) = o'mi' ^'^(r)^'^ {<!>)• Appendix 
IeI we give the explicit forms for the shear tensors a^i, and the final form of a(r)(<i>)- These calculations are more lengthy than 
the same calculations using the Boyer-Lindquist coordinate but straightforward. Here, we simply express the shear rate <y(r)(4,) 
as 
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du^ dl , , 

cr(,)(^) =a = cr, + a„— +cr, — . (64) 



From Eqs. (|60|l . (|63|l and (|64p . we can derive the equation for dl/dr having the singular point which we call viscous point. 
Based on the shear stress and the shear rate calculated in the last section, the dissipation function $ is calculated as 

$ = -a^^t"" = -cr(c,)(/3)t<"'(''' = -2cr(^)(<^)t('''('*', and finally we obtain 

$ = -2oS, (65) 
where a and S are given by Eqs. ((60} and (|64|l . 



4 BOUNDARY CONDITIONS 

The accretion flows plunging into the black hole supersonically must pass the sonic point where the accretion velocity become 
larger than the sound speed. On the other hand, when the causal viscosity prescription is used, the accretion fiows pass the 
viscous point where the accretion velocity become larger than the speed of the viscous diffusion. In order to smoothly pass 
the sonic point and the viscous point, the fiows must satisfy the boundary conditions at the sonic point ( H4.f |l and the viscous 
point ( H4.2p . By using the causal viscous prescription, the boundary conditions at the outer regions or the inner regions of 
the accretion flows are not required. The boundary conditions at the event horizon which are used in some past studies are 
not required in the present study. 



4.1 Boundary Conditions at the Sonic Point 

In order to obtain the boundary conditions at the sonic point, we need the equation which do not contain the derivatives 
except du^ /dr. 

The pressure p and the internal energy u of the accreting fluid is usually a function of the rest-mass density po and/or 
the temperature T, i.e. p — p{po,T) and u = u{po,T). Then, the derivative dp/dr is calculated as 

dp ^ / dp\ dpo ( dp\ dT 

dr \dpoJ^ dr ^\dT),„ dr ' ^''''> 

dr ^ \dpo )^ dr ^ Urjp„ dr- ^ ' 

In the case of the pressure and the internal energy is a function of the rest-mass density only, all the thermodynamic quantities 
can be written by the rest-mass density only, i.e. {dp/dT)pg = {du/dT)pg = 0. In such case, with Eqs. (|66p . (|67|l and the mass 
conservation given by Eq. ((TSJ, we can obtain the derivative dp/dr described by the derivatives du^ /dr and dl/dr as 

dlnp „ „ du'' „ dl , . 

dr dr dr 

Usually, the derivation for the boundary condition at the sonic point become lengthy and complex for the general relativistic 
accretion flows. Eq. (|68|l is the key equation in order to simply treat the boundary condition at the sonic point. In the case 
of the general equation of state where the pressure is a function of both the rest-mass density and the temperature, if we use 
the energy equation given by Eq. (I52p . we can also have the derivative dp/dr with the form as Eq. (|68|l . The examples of Pr, 
Pu and Pi for the isothermal disk, the polytropic disk and the ADAF with the general relativistic equation of state are given 
in Sec. [7] From Eqs. (|32|) . (|64|l and (|65|l . the term nm containing the effects of the heat inertia is calculated as 

TlHI = nHI + TlHI + TlHI ^ , (69) 

where 

T ^ U U S £ U S /r^r\\ 

riHi = o",,, nHI = fu, ^Hi = <Ti- (70) 

pan poll poV 

Here, 5* is calculated from Eqs. (|38p and (|58|) as S* = —pou^lrjl — j)/ F. If the heat inertia effects by the flux is required, we 
change ri^i and jihi according to Eq. (|3I|) . On the other hand, T)^, and AC are defined in the next section in order to 
have dl/dr with the form described as 

dl dv^ 

v^^=n:+n:^. (71) 

dr dr 

By substituting dp/dr and nni described by Eq. (|68|l and (|69[) . respectively, into the radial momentum equation given by 
Eq. (|33}, the derivative du^ /dr can be calculated as 
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where 

Vs = u^^h^^cl(^P^ + P,^'^-(n^, + ni,,^y (73) 

J^s = nacc + fe^'-c?(p. + ft^) + (n^i + n^Hi^) . (74) 

Here, we use Eq. (|7ip to remove the derivative dl/dr. 

In order to pass the sonic point smoothly where = 0, the condition Na ~ must be satisfied at the sonic point. So, 
the boundary conditions at the sonic point are 

v,=^^s = 0. (75) 



4.2 Another Boundary Condition 

4-2.1 For Type A Causal Viscosity: Boundary Condition at Horizon 

For the type 1 causal viscosity prescription, we put the boundary condition at the horizon in order to vanish the angular 
momentum transportation at the horizon. At r = r+, from the condition v = 0, rji = j is required and is used as the boundary 
condition. In this case, the parameters given at the sonic point determine the transonic solution. The difi'erential equation 
for the angular momentum, i, is calculated from the equation for the angular momentum conservation and the Navier-Stokes 
prescription of the viscosity as t^„ — —2pQrji/a^,^ as 



where 








r 

U ;^ 


^ / tr 1 rr , r4> \ r 

= 2^9 9t<p,r + g gr4,,r + g "^gu.rju ~ 




(77) 


U4,-t 


1 




(78) 




1 




(79) 




= ^{gt.j>,r + g,i„j,,r^)u , 




(80) 


and Qr 


is given in Appendix. Here, the shear rate 


cr'"^ is calculated as 




_ u'~{ril- j) 




(81) 


^ 







when the value of the kinematic viscosity coefficient is not zero, i.e., v ^ 0. When v = 0, this shear rate is zero, i.e. a^^ = 0. 



4-2.2 For Type B Causal Viscosity: Boundary Conditions at Viscous Point 



The boundary conditions at the viscous point are calculated from the equation for dl/dr. In the same way as Eq. (I68|l . the 
derivatives of rj, Hg and F with respect to r are described by the combinations of du^ /dr and dl/dr as 



d In 77 
dr 
dlnHg 
dr 
d\nF 
dr 



Hr + Hu 



dl 



du/_ 

du'' 



Hi 



dr dr 

du'' d£ 
Fr + Fu—+Fe- 



dr dr 

By substituting Eqs. ([Ml), IMl), (US), dM]) and (dH) into Eq. (O and using the relation 5" 
equation for dl/dr containing no derivatives except du^ /dr can be calculated as, 

dI_N^ 
dr T>v ' 

where JVy = TV^'' + JV^ (du'' / dr) and 



-pou''{ril - 



r\2 



K) 



2aiFcl 



+ 



2(7/ Fc: 



l+lr)e- 
Irjr- 



Qe 



SF 
SF 

Pom'' 



Gr 



(82) 
(83) 
(84) 

Qi)/F, the 
(85) 

(86) 
(87) 



12 Rohta Takahashi 



AC 



0"u 



Qu ^ SF 
i r<J"u 



2atFcl 

Here, Gr, Gu and Ge are defined as 

G, = 2/r + + Gu = Fu+ H^, Ge = Fe 

in order to iiave tlie relation 

M 



2 dlnF dlnHg _ du" 
h ; \ ; — Gr + Gu — ; — 

r ar ar ar 



Gi 



dr ' 



in the denominator of Eq. (|60|l . and Qr, Qu and Qe are determined in order to iiave the relation 

dQi _^ , A du^ ~ dl 
dr dr dr 



(89) 



(90) 



(91) 



In this study, we neglect the angular momentum loss by the radiation, i.e. we set Qr = Qu = Qe = 0. 

In order to pass the viscous point smoothly where I?„ = 0, the condition Af^ = must be satisfied at the viscous point. 
So, the boundary conditions at the viscous point are 



Vu=J^v= 0. 



(92) 



5 COUPLED DIFFERENTIAL EQUATIONS TO BE SOLVED 

For the general equation of state, the transonic solutions are obtained by numerically solving the coupled differential equations 
for the dynamic variables, e.g. u*", £, and the thermodynamic variables, e.g. T. In the case of the special thermodynamic 
relations, such as the isothermal flows and the polytropic flows, the thermodynamic variables can be calculated from the 
dynamical variables. In these cases, we only solve the coupled differential equation for the dynamic variables. In this study, 
we treat the radial component of the four velocity, u' , and the angular momentum, £, as the basic dynamical variables to be 
solved. That is, for the case of the special thermodynamic relations where the thermodynamic variables are calculated from 
the dynamic variables, we solve the coupled differential equations for and £ described as 



dr Vs ' 

d£ Nu ( AC dit^AC 



dr Dv \ "Dv dr "Dv 



(93) 
(94) 



In the following sections where the transonic solutions for the isothermal flows and the polytropic flows are calculated, we 
solve these two differential equations. 

On the other hand, for the general equation of state, the differential equation for the thermodynamic variables is usually 
solved in addition to Eqs. (|93p and (|94|l . The differential equations for the thermodynamic variables are derived by using the 
energy equation given by Eq. (I52p . In this study, we treat the temperature T as the basic thermodynamic variable whose 
differential equation is numerically solved. The other thermodynamic variables, such as the rest-mass density po, the sound 
velocity Ca, are calculated from , I and T by using the mass conservation equation given by Eq. (|18p and the equation of 
state. Here, we derive the general form of the differential equation for T by using the energy equation. By differentiating the 
mass conservation given Eq. (|18|) . the disk thickness Hg — Csr/£, and the sound velocity Cs — [p/('?Po)]^''^ with respect to r, 
we obtain 

2 ^ dlnpo ^ rflnj^n _^ dlnHg ^ ^ 
r dr dr dr ' 

dlnHg dine. 1 1 



dr dr r £, dr ^ 

^dlncs _ dlnp dlnpo ^ dln-q 
dr dr dr dr 

Here d£,/dr can be written by the linear combination of du^ /dr and dl/dr as 

d£* „r „7, du^ „e d£ , . 

dr dr dr 

where the coefficients £1, and di can be calculated analytically or numerically. We newly define 

-.= fel' "'fe)/ 
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Here, rjp and r]T are related to Pp, Pt, Up and Ut as rjp = {uUp + pPp — u — p)/{rjpo) and rir — {uUt + pPr) / {ripo). From 
Eqs. ((66j, (|67t . (|95)) . ((96)) and ((97)) . we obtain the differential equation for T and po as 

^ ^ T.+T„^+T.J, (100) 
dlnpo du^ dl 

Here Tfc and pk (k — r, u and ^) are calculated as 

Tfe = CfA'fe+Cjg^K, (102) 
Pk = CfA-fc+Ca^q^/ii", (103) 
where Xk {k = r, u and £) are given as 

Q 1 /" 

and Afy, A'y, A'p and A'p are given as 

2 / u + p\ T _}_tp nP 2f/T p {Pt - Vt) , . 

Here, Xd is defined as 

A-c = (Pp - r,p + 1)Ut - {Pt - r,T) (Up - (106) 

On the other hand, {k — r, u and £) are defined to satisfy the relation 

+ - ± ±dvr -j-d£ 

By using the coefficients Tk and pfe {k = r, u and i?) calculated above, the coefficients Pk, rjk and Hk {k — r, u and £) are 
calculated as 

Pk = -PTTk-PpPk, (108) 

= VTTk + VpPk- (109) 



From Eq. (|95}, we obtain 

Hr^-^-P^, H^^-\-^, m = (110) 

r- po M po po 

For the general equation of state, the transonic solutions are obtained by solving the differential equations of , I and T 
given by Eqs. (|93p . (|94|) and (|100|) . In the later sections, for the ADAF with the general relativistic equation of state and the 
supercritical accretion disk, we solve the differential equations for , I and T. Since the rest-mass density po is calculated 
from the mass conservation equation given by Eq. (|18() . we do not solve the differential equation for po. 



6 CALCULATION METHOD 

By using the formalism developed until the last sections, we solve the coupled differential equations to obtain the transonic 
solutions. The calculation method for the transonic solutions are not unique, and actually, past studies use several method. 
Here, we show one of the calculation methods to obtain the transonic solutions. 

(i) First, we tentatively choose some value of (or as^s) for given values of and j, and calculate u^, £s, {du^ /dr)s, 
{dl/dr)s and {dT/dr)^. Here, the differential values at the sonic point are calculated by using the L'Hopital's rule. 

(ii) Next, we solve the solutions in the range < r < r„. In order to do this, we solve the coupled differential equations for 
, £ and T from the sonic point to the viscous point by using, e.g., the Runge-Kutta algorithm. Usually, for the initially selected 

value of Ts (or Us.s), the calculated solution does not pass the viscous point where two boundary conditions D„ = A/'i, — 
are satisfied. In such case, we return to step 1 and again choose the different values of (or as,s) for given values of and 
j. After repeating these procedures, we can determine the value Ts (or as,s) which gives the solution satisfying the boundary 
conditions Ds ~ A/'s = at r = rs and Vv — Nv = Q aX r = Tv 

(iii) After solving the solutions in rs < r < ri,, we solve the coupled differential equations in the range r„ < r by using the 
values of Ts (or as,s) for given values of rs and j by using, e.g., the Runge-Kutta algorithm. 

(iv) Finally, we solve the coupled differential equations in the range r < rs by using the values of Ts (or as,s) for given 
values of Vs and j by using, e.g., the Runge-Kutta algorithm. If the solutions are connected with the horizon as usual solutions. 
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we can solve the solutions inside the horizon. On the other hand, if the solutions are not connected with the horizon as the 
alpha-type solutions, the numerical integrations are stopped before the horizon because there is no stationary solutions for 
such parameters of and j. 

The third step and the fourth step can be interchanged. By this procedure, the transonic solutions are obtained for given 
values of ra and j without the boundary conditions for the inner regions (r < Ts) of the outer regions (r„ < r). By using these 
procedures, we can basically cover all the possible values of and j. That is, by these methods, in principle, all the possible 
stationary transonic solutions can be calculated because we can use all the possible sonic point, and the transonic solution is 
calculated from the sonic point. 



7 APPLICATIONS AND SAMPLE SOLUTIONS 

In this section, we give the numerical solutions for the ideal isothermal accretion flow f ^7.1|l . the polytropic disks ( ^7.2[l . the 
ADAF with relativistic equation of state ( ^7.3|l and the adiabatic accretion flow (i l7.4p and the formulation for the supercritical 
accretion flow fi 37.5p . Based on the thermodynamic relations, we first calculate the coefficients Pk, rjt and Hk {k = r, u and 
£) which are used in the calculations of du^ /dr and dl/dr. In addition, we calculate {k = r, u and t) if required. For the 
ideal flows, we only solve the differential equation of . For the viscous polytropic flows, we solve the coupled differential 
equations of and I. For ADAFs with the relativistic equation of state, in addition to the differential equation of vT and I, 
we also solve the differential equation of T simultaneously, which is derived by using the energy equation. In this section, for 
the viscous solutions, the kinematic viscosity v is assumed to described by the alpha viscosity as v = a„Cs/n. 



7.1 Application 1 : Ideal Isothermal Accretion Flow 

By using the formalism developed until the last sections, we first shows the numerical solutions of the horizon-penetrating 
solutions for the ideal isothermal accretion flow which is one of the simplest transonic accretion flow. Here, we only solve the 
differential equation for it"" by assuming constant specific angular momentum I and sound speed Cs, and 77 = 1 is also assumed. 
We use the coefficients Pr, Pu and Pe described as 

Pr=-, Pu = ^, Pl=0. (fff) 

r 

By substituting Eqs. into Eqs. ((72)), the differential equation for for the ideal isothermal flows are obtained. We 

numerically solve this differential equation for and obtain the transonic solutions. When calculating the numerical solutions, 
the rest-mass density po is determined from , I and Cs by using the mass conservation equation given by Eq. psp . For the 
ideal isothermal flows which are solved in this section, since we assume constant I, there is no viscous point in the global 
solution of the transonic accretion flow. 

In Fig.[T] in the parameter spaces rs-c^ we plot lines of constant critical values of A(= ~u^/ut). The critical sound speed 
is plotted for non-rotating (a/m — Odeft panel) and rotating (0.95: right panel) black holes. We calculate transonic solutions 
with the critical values plotted by the filled triangles in Fig. [T] The resultant transonic solutions for non-rotating black holes 
are plotted in Fig. [5] We also give the transonic solutions for rotating black holes in Fig. O For both Fig. [5] and Fig. [31 
the transonic solutions calculated by using the Kerr-Schild coordinate {left column) or the Boyer-Lindquist coordinate [right 
column) are plotted. The angular velocity Q, and directly reflect the effects of the coordinate singularity when we use the 
Boyer-Lindquist coordinate. That is, for the solutions calculated in the Boyer-Lindquist coordinate, the angular velocity, Q,, 
is equal to the angular velocity of the frame dragging, lu, at the horizon, and is diverged at the horizon. These features are 
clearly seen in both Fig. [5] and Fig. O We also show Q, ~ u in the inserted box in the panel showing in the right column 
of Fig. 13] The feature that the angular velocity of the accretion flow written by the Boyer-Lindquis t coordinate is equal to 



the angular velocity of the black hole lo = a/2mr at the horizon is pointed out by iKomissarovl 1 20041 ) who also found that it* 
remains finite at the event horizon and differs from the angular velocity of the black hole in Kerr-Schild coordinate. The 
coordinate singularity in Boyer-Lindquist coordinate is also relevant to the feature that the world lines of Boyer-Lindquist 
LNRF become null on the event horizon and thus cannot correspond to any physical observer. For the outside region of the 
horizon, the lines for u*" are same for both calculations using Kerr-Schild coordinate and the Boyer-Lindquist coordinate. We 
plot two types of transonic solutions which have the sonic radius in the inside region or the outside region. These two types 
of solutions correspond to the solutions named type I and type II in Peitz & Appl (1997). Similar solution patterns are also 



obtained bv lFukud (|l987|) 



For the accretion flows calculated by using the Boyer-Lindquist coordinate, we also plot the results for the flows which 
are firstly calculated based on the Boyer-Lindquist coordinate and then transformed to the flows written by the Kerr-Schild 
coordinate by the transformations of four velocity given by Eqs. (|10|l and . These results are plotted by the short dashed 
lines in the right panels for Q, and it* of Figs. [2] and Fig. [3] These solutions outside the event horizon are same as those 
calculated based on the Kerr-Schild coordinate presented in the left panels of Figs. [2] and Fig. [S] 



Horizon-Penetrating Transonic Accretion Disks around Rotating Black Holes 15 




Figure 1. Parameter spaces rg-c^ showing lines of constant critical values of A(= —u^/ut). The critical sound speed is plotted 
for non-rotating {a/m = O.left panel) and rotating (0.95: right panel) black holes. For non-rotating black holes, contours correspond to 
A = 2.0, 3.2, 3.4, 3.6 and 4.0 (from right to left). For rotating black holes, contours correspond to A = 1.5, 1.9, 2.1, 2.2, 2.3 and 2.5 (from 
right to left) . Transonic solutions with the critical values plotted by the filled triangles are calculated in Fig. [2] for non-rotating black 
holes and in Fig. [3] for rotating black holes. 



Outside the horizon, from the results for the accretion flow calculated by using the Boyer-Lindquist coordinate shown in 
the right panels of Figs. [2] and [S] we can obtain the solutions given in the left panels for the Kerr-Schild coordinate by using 
the transformation given by Eqs. (|10p and (|lip . While for the ideal accretion flows the accretion flows calculated by these two 
procedures have same results, for the viscous flows the solutions calculated by these two procedures do not have the exactly 
same results. See the discussion in the last section. 



7.2 Application 2 : Polytropic Accretion Flow 

Here, we show the transonic accretion flows with the polytropic equation flows. Although this equation determine the general 
energy equation, in this paper, we only consider the accretion flows with the polytropic equation of state as 

p = Kpl, (112) 

where K is constant and F is the adiabatic index (or the ratio of speciflc heat). The adiabatic index F is related to the 
polytropic index A'' as F = 1 -I- 1/N. The internal energy u is given hy u — p/{r — 1). When solving the transonic solutions for 
polytropic accretion flows, we do not use the energy equation in Eq. (|52p . All the thermodynamic variables are expressed by 
the rest-mass density po for given adiabatic index F and the constant K. The relativistic enthalpy rj [= 1 + {u + p)/po] and 
the sound speed [= p/ipov)] are calculated as 

(114) 

(115) 
(116) 
(117) 



P 



po + [F/(F-l)lp- 

From Eqs. (|112p . (|113p and (|114p . dp/dr, drj/dr and dcs/dr and can be calculated from dpo/dr as 

dp Tp dpo 

dr Po dr ' 

drj Tp dpo 

dr Pq dr ' 

dcs _ (F — l)c3 dpo 
dr 2r]po dr 

From the disk thickness given by Eq. (|51|l and the mass conservation given by Eq. (|18p . we obtain Eqs. (|95|l and (|96fl . By 
eliminating dHe/dr from Eqs. (|95|l and (|96p and substituting Eq. (|117p . the derivative dpo/dr is calculated as, 

dpo -2r]po 





"4 £1 


/I C\ 


du'' 


£i dt 


1 


r^ T,^ 




dr 


~ J^dr. 



dr 2-n + r - 

By substituting Eq. pisp into Eq. pisp . Pr, Pu and Pe are calculated as, 



(118) 



Pr ^1 



/4 r\ / 1 £ 
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Figure 2. Sample transonic solutions for the ideal isothermal flows when a/m = 0. We plot transonic solutions with the critical values 
A = 2.0, 3.2, 3.4, 3.6 and 4.0 (from right to left) with the sound speed = 0.03. The transonic solutions calculated by using the 
Kerr-Schild coordinate {left column) or the Boyer-Lindquist coordinate (right column) are plotted. The radius of the horizon is denoted 
by the dashed line, and the angular velocity of the frame dragging is plotted by the dotted lines in the panel showing f2. 
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Figure 3. Sample transonic solutions for the isothermal ideal flows when a/m = 0.95. We plot transonic solutions with the critical 
values A = 1.5, 1.9, 2.1, 2.2, 2.3 and 2.5 (from right to left) with the sound speed = 0.1. The transonic solutions calculated by using the 
Kerr-Schild coordinate [left column) or the Boyer-Lindquist coordinate {right column) are plotted. The radius of the horizon is denoted 
by the dashed line, and the angular velocity of the frame dragging is plotted by the dotted lines in the panel showing Q. We also show 
n — oj in the inserted box in the panel showing O in the right column. 
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Figure 4. Sample transonic solutions for tlie polytropic ideal flows for a/m = 0. We plot transonic solutions with the critical values 
A = 2.0, 2.6, 3.0, 3.2, 3.4, 3.6 and 4.0 (from right to left) with the critical sound speed = 0.05. The radius of the horizon is denoted 
by the dashed line, and the sonic points are plotted by the flUed triangles. 




Here, we define Xp = 2riT/{2r] + F — 1). In the same way, rjr, rju and rje are calculated as, 

Vr = -CsPr, rju = -c^Pu, Tje = -c^Pi, (120) 
and Hr, and Hi are calculated as, 

By substituting Eqs. (|119|l . (|120|l and ^211) into Eqs. (f72)) and the differential equations for u"" and £ for the 

polytropic accretion flows are described as du^ /dr — Ns/'Da, dljdr = Ml jT)^ + (Ml /'D^)dvJ' / dr which are the basic coupled 
differential equations to be solved. In this section, we calculate the ideal polytropic flows and the viscous polytropic flows. For 
the ideal polytropic flows, we assume "ql =constant and only solve the differential equation for , and the global transonic 
solutions do not have the viscous point. For the viscous polytropic flows, we solve the coupled differential equations of u"" 
and I, and obtain the transonic solutions satisfying the boundary conditions at the sonic point and the viscous point. For all 
the viscous flows, we assume the alpha viscosity — 0.1. The radius of the sonic point is determined in order to pass the 
viscous point. In the same way as the previous section, when calculating the numerical solutions, the rest-mass density po is 
determined from , £ and Cs by using the mass conservation equation given by Eq. (|18p . 

We show the sample solutions for the horizon-penetrating transonic solutions of the ideal polytropic flows for a/m — 
and 0.95 in Fig. |4]and Fig. [5l respectively. The sonic points are plotted by the filled triangles, and the radius of the horizon 
is plotted by the dashed lines. The transonic solutions for the viscous polytropic flows for a/m = and 0.95 are given in Fig. 
[6] and Fig. [T] respectively. For the viscous flows, the sonic points and the viscous points are plotted by the filled triangles and 
squares, respectively. All the viscous polytropic solutions presented here become super-Keplerian flows in the outer region, 
and the sound speed diverged. The outer region of these solutions correspond to the thick disk solutions. 



7.3 Application 3 : ADAFs with Relativistic Equation of State 

Here, we calculate the transonic solutions for the ADAF with relativistic equation of state, where the energy equation given 
by Eq. ((52} is required in order to close the coupled differential equations. For the Boyer-Lindquist coordinate which have the 
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Figure 6. Sample transonic solutions for the viscous polytropic flows when a/m = 0. We plot transonic solutions with the constant 
specific angular momentum j = 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5 and 3.7 (from right to left) with the critical sound speed = 0.05. The 
radius of the horizon is denoted by the dashed line. The sonic points and the viscous points are plotted by the filled triangles and squares, 
respectively. 



coordinate singularity at the horizon. iGammie fc Popham l|l998t ) and lPopham fc Gammiel (|l998l ) solve the transonic solutions 

by using the causal viscosity prescription. 

Here, we use the equation of state same as lGammie fc Popham (Il998l). The pressure p and the internal energy u is given 
by the rest-mass density po and the temperature T as (jChandrasekhail 19391 : ICox fc Giulilll968h 



u = po 



3T + 



Ki{l/T) 
K2{l/T) 



(122) 
(123) 



wh ere Kn 's are the modified B essel functions of the second kind of order n. The internal energy u is well fitted by the function 
as ijGammie fc Popham|[l99i ) u = poTg(T) where g{T) = (45r^ + 45T + 12)/(15r^ + 20T + 8). The relativistic enthalpy r) 
and the sound velocity Cs become a function of the temperature T as, 

0. T 



, = 1 + T[1 + ,(T)], cf = r+T[i + ,(T)]- 
Here, rj and Cs are functions of temperature, T. By using these equations, we can obtain 



, d\ng(T) 



1 / dlng(T)\ 

1 + U + u 

Po?7 V dluT ) 



d\nT 

On the other hand, by using the assumption q^^^ — 0, we can calculate the coefficients (k = r, u and £) as 



1^ 



(124) 



(125) 



(126) 
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By substituting equations given by Eq. (|125p into Eqs. (|102|l . (|108|) and (|110p . we can calculate the coefficients Pk, Hk, rik 
and Tfc {k — r, u and £). Now, it is noted that the relativistic enthalpy rj is also calculated by the equation in Eq. (I124|l . Then, 
the derivatives du^ /dr, dl/dr and dT/dr are obtained and numerically solved in order to calculate the transonic solutions for 
the adiabatic accretion disks. In Fig. [S] we show the sample numerical transonic solutions for the ADAF with the relativistic 
equation of state. Sample numerical solutions are calculated for a/m — 0.0, 0.5, 0.95 and 0.99999 with Ts — 0.1, and we plot 
the four- velocity components, and , the angular velocity fi, the dimensionless temperature T, the relativistic enthalpy ri 
and rj£. The solutions are calculated so as to satisfy zero shear stress at the horizon. The positions of the horizons and the 
sonic points are denoted by the blank circles and the filled triangles, respectively. The dashed lines in the panel for rji {bottom 
right) are the angular momentum for the Keplerian motion when a/m = 0.0, 0.5, 0.95 and 0.99999 (right to left). All solutions 
pass the event horizon smoothly and have nearly same fiow patterns in the outer region. As denoted in the previous section, 
the accretion flow is plunging into the black hole with the angular velocity which is different from the angular velocity of the 
black hole's rotation at the event horizon. While the accretion fiow have the smaller values of lu""! for larger values of black 
hole spins, the gamma factor 7 = au* become larger for the larger values of the black hole spins. This is because of the effects 
of the black hole rotation enhances the rotational velocity of the accretion flow, i.e., for larger black hole spins the accretion 
flows have the larger values of the angular three velocity. The relativistic enthalpy rj = 1 + (u + p) /po of the hot accretion flow 
is in general larger than unity near the horizon. In the sample solutions for a/m = 0.95 and 0.99999, the accretion flows are 
sub-Keplerian in all the region. On the ot her hand, for a/m = 0.5, the middle part of the accretion flow is super-Keplerian. 
The similar feature is also pointed out bv lPeitz fc Appll l| 19971 ). 
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Figure 8. Sample numerical solutions for transonic accretion flows of ADAFs with relativistic equation of state when a/m = 0.0, 0.5, 
0.95 and 0.99999 with Ts = 0.1. We plot the four-velocity components, •u'" and «*, the angular velocity Q, the dimensionless temperature 
T, the relativistic enthalpy rj and rji. The solutions are calculated so as to satisfy zero shear stress at the horizon. The positions of the 
horizons and the sonic points are denoted by the blank circles and the filled triangles, respectively. The dashed lines in the panel for 
ri£ {bottom right) are the angular momentum for the Keplerian motion when a/m = 0.0, 0.5, 0.95 and 0.99999 (right to left). Here, we 
assume the alpha viscosity parameter 0.01 and type A causal viscosity. 
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7.4 Application 4 : Adiabatic Accretion Disk and Standard Accretion Disk 



In this section, we calculate the accretion flows where the viscous heating rate is balanced with the radiative cooling rate, 
Ivis ~ ^rad- Thls assumptiou is usually used in the calculations of the standard accretion disks. From this assumption and the 
energy equation given by Eq. ([55}, we can show that the entropy change of the accretion flow is zero, i.e. dS = 0. So, here, 
we call the disk with this assumption as adiabatic accretion disk. 

Here, we calculate the transonic solutions for the adiabatic accretion disks whose mass accretion rate is near or beyond 
the Eddington mass accretion rate. For the supercritical accretion flow the contribution from the radiation pressure of photons 
can not be neglected. When the specific heat at the constant volume, cv, is inde pendent of the tempe rature, the pressure p 
and the internal energy u for the flows containing gas and radiation are given as I Chandrasekhai 193y), 



P 



Pa +Pr, 

1 

Pa 



7-1 
where pg and 



Pa = 



PoT, 



3pr, 

are the gas pressure and the radiation pressure, respectively, described as 

1 ^4 
Pr = garadJ • 



(127) 
(128) 

(129) 



Here, 7 is the ratio of the specific heats, fcs is the Boltzmann constant, /i is the mean molecular weight, ruH is the Hydrogen 
mass, flrad is the radiation constant, and we use the dimensional representation of the pressures. These assumptions are 
sometimes used in the past studies for the standard accretion disk (e.g. Shakura & Sunyaev 1973, Novikov & Thorne 1974, 
Page & Thorne 1974, Matsumoto et al. 1984). For fiows with dS = 0, the energy equation, u'' poT{dS/dr) = g+j^ — q^^^{= 0), 
can be given by using the generalized adiabatic exponents as 



Fa - 1 
where 



dp _ P_rfpo 



dr 



po dr 



i(= 0), 



Fi 
F3 



(4-3/3)2(7-1) 



1 + 



/3 + 12(7- 
Fi -/3 
4-3/3' 



(130) 

(131) 
(132) 



Here, /3 is the ratio of the gas pressure to the total pressure, i.e. /3 = Pg/p- In the calculations of the transonic fiows, we use 
the dimensionless pressure, internal energy, rest-mass density and the temperature. By using the dimensionless prescription, 
the pressure and the internal energy can be calculated a,s p — pg + pr and u — Pg/ij ~ 1) + Spr where pg = poT and pr — r*/3 
where the temperature is normalized by mp(? /ks, the rest-mass density is firstly normalized by {ar^d/ c?){mpC? /ks)*' to the 
dimensionless quantities and is secondly normalized so as to satisfy M = 1. Here, is the proton mass and c is the speed 
of light. Since the left-hand-side of Eq. (|130p derived from the condition dS = is equivalent to that of Eq. (|55p . we have the 
relations 



\dTjf 

du 
dpo 



1 



1 



1 \dT), 

dp 
dpo 



F3-I 

From the equation of state and Eqs 
Pp = /3, Pt =4-3/3, Up 



-Fi- 



po 

lfT33l) and 



po 

(fT34l) 



we have 



2(4 - 3/3) + 



/3 



7-1 



Ut 



P / 4-3/3 \ 
u \r3-l) 



(133) 
(134) 

(135) 



where u/p = /3/(7 — 1) + 3(1 — /3). By substituting equations in Eq. (|135p into Eqs. (|102|) . (|108p and (|110|) . we can calculate 
the coefficients Pk, Hk, rjk and Tk {k = r, u and £). Then, the derivatives du^ /dr, dl/dr and dT/dr are obtained and 
numerically solved in order to calculate the transonic solutions for the adiabatic accretion disks. In Fig. |9l we show the 
sample transonic solutions of the adiabatic accretion disks with the black hole mass A^bh = IOM0 and the mass accretion 
rate M = lOA/Edd where A^Edd = 1-4 x 1O^^(A/bh/A/0) [g s~^] is the Eddington mass accretion rate. Sample solutions are 
calculated for a/m = 0.0, 0.5, 0.95 and 0.998 with dimensionless temperature Ts — 0.05, 0.01, 0.05 and 0.2, respectively. In 
Fig- El the radial component of four velocity (left panel) and the dimensional temperature T [K] {right panel) are plotted. 
The positions of the horizons and the sonic points are denoted by the blank circles and filled triangles, respectively. The 
solutions are calculated so as to satisfy zero shear stress at the horizon. The temperatures of all solutions become smaller by 
several order inside the marginally stable orbit. This feature is same as the standard accretion disks. On the other hand, inside 
the marginally stable orbit, the absolute value of the radial component of the four velocity become large. This is because the 
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Figure 9. Sample numerical solutions for the adiabatic accretion disks when a/m = 0.0, 0.5, 0.95 and 0.998 with dimensionless 
temperature Ts = 0.05, 0.01, 0.05 and 0.2, respectively. The black hole mass Mbh = IOMq and the mass accretion rate M = WM^^^ 
are assumed. The radial component of four velocity and the dimensional temperature T [K] are plotted. The positions of the horizons 
and the sonic points are denoted by the blank circles and filled triangles, respectively. The solutions are calculated so as to satisfy zero 
shear stress at the horizon. Here, the alpha viscosity parameter 0.01 and type A viscosity are assumed. 

accretion flow approach the free fall motion with angular momentum in this region and plunge into the black hole horizon 
with high value of the gamma factor. 



7.5 Application 5: Supercritical Accretion Disk with Photon-trapping Effects 

Based on the formulation described above, we also calculate the simplest version of the transonic solutions of the slim disk 
of the supercritical accretion flow where the mass accretion rate is larger than the super-Eddington mass accretion rate. 
In general, for the supercritical accretion disks, the assumption of adiabatic changes, dS — 0, is not valid because in the 
vicinity of a black hole the pressure gradient enhance the radial velocity of the flow and the advection term in the energy 
equation can not be neglected. In such case, the energy equation with the effects of advection cooling, the radiation cooling 
and the viscous heating should be solved. In addition, near the horizon, photons are trapped within matters and can not 
escape from the accretion flow, and the flow become advection dominated state. The past studies actually solve the transonic 
solutions of the supercritical accretion flow (e.g. Watarai et al. 2001, Watarai & Mineshige 2001, Shimura & Manmoto 2003) 
by assuming that the speciflc heat at the cons tant volume, cy, is independent of the temperature. For supercritical accretion 
flow, the heat inertia can not be neglected ( Beloborodov 19981 ) and the photon-trapping effects near the horizon is also 



important. For the Boyer-Li ndquist coordinate which have the coordinate singularity at the horizon, Beloborodov ( 19981 ) and 



Shimura fc Manmotol (|2003l ) solve the transonic solutions based on the acausal viscosity prescription 



From here, we just show the formulation of the transonic solutions for the supercritical accretion flows with effects of 
the heat inertia and the photon trapping with the general form of the specific heat at the constant volume, i.e., here, we do 
not assume that the specific heat at the constant volume, cv , is independent of the temperature. For such fiows, the internal 
energy u is calculated as 

u=p,g{T) + 3pr, (136) 

where g{T) is the same function used in the calculations of the transonic solutions of ADAFs with relativistic equation of state 
described in the previous section. First, we roughly estimate the photon trapping effects around black holes. In the optically 
thick region for photons in the supercritical accretion flows, the photons can escape from the disk surface after the diffusive 
processes in the disk. This holds only when the radiative diffusion timescale is shorter than the accretion timescale (see e.g. 
Katz 1977, Begelman 1978, Ohsuga et al. 2002, Kawaguchi 2003). The diffusion velocity is roughly calculated as «diff ~ c/(3r). 
Then, the diffusion timescale of photons produced at equatorial plane is written as t^ig — H/v^iff. On the other hand, the 
accretion time scale is tacc = r/{—u^). When tdis > tacc, photons are trapped in the accretion disk and plunged into black 
hole without escaping from the disk surface. From the condition tdut > tacc, we can derive the photon-trapping radius, r'^ap, 
within which the parts of photons begin to be trapped as rj^ap ~ [3R/{2TTc)]HeM where R is the mean opacity for photons. 
The radiation term including effects of the photon-trapping, the electron scattering and the free-free absorption is calculated 
l7ad ~ /trapQ^d g whcre q~^^ g IS the cooling rates when no effects of photon trapping. Now, /trap ~ 1 means no photon 
trapping and /trap = means complete photon trapping. In this study, we assume /trap = 1 for r > rtrap, but /trap ~ 
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Figure 10. Transonic solutions of ADAF calculated by using the Kerr-Schild coordinate {silid lines) and the Boyer-Lindquist coordinate 
(dashed lines) for a/m = 0.95 and 0.99999. The parameters are same as Fig.|8] The blank circles show the position of the horizon. 



for r ^ rtrap, for simplicity. The radiative cooling term, q~^^ ^ [erg cm""^ without the photon trapping is calculated as 

9rad ~ /{2rHg) where F~ [erg cm~'^ s~^] is the energy loss rate from the disk surface which is calculated by using the 
Rosseland ap proximation as F~ = [IQassT'^) / (SKporHg) in the outside region, and F~ is calculated by the free-free emission 
inside region ( Beloborodov 1998h . Here, (Jsb is the Stephan-Boltzmann constant. We can calculate the coefficients (k = r, 
u and £) as 



-ISGt 



9radi 9u 



-2Sau, = —2Sai 



On the other hand, from the EOS given by Eqs. (|127|l and (I128p . we obtain 



Pp = P, Pt=4 



3/3, Up = 



+ 12(1-/3) 



(137) 



(138) 



where u/p = /3p(r) + 3(1 - /3) and dg/dt = 15(15r^ + 24T + 8)/(15T^ + 20T + 8f. By substituting Eqs. (IT38|l and l{T37\ 
into Eqs. (|102p . (I108|l and (1110(1 . we can calculate the coefficients Pt, Hk, rjk and Tk (k — r, u and £). Then, the derivatives 
du^ /dr, dl/dr and dT/dr are obtained and numerically solved in order to calculate the transonic solutions for the supercritical 
accretion disks with effects of advection cooling. 



8 CONCLUDING REMARKS 

Before summing up the results of this study, it may be better to note the arguments about the causality of the viscous flows at 
the event horizon. As already pointed out by Popham & Gammie (1998), the formalism of the causal viscosity prescription give 
the solutions with the finite values of the outward energy flux and the outward angular momentum at the horizon, and they 
stated that this property does not suggest the causality violation (Popham & Gammie 1998) . If the viscosity is described by the 
fluctuations of the Maxwell stress and/or the Reynolds stress from the mean values such as the angular momentum transport 
by the magnetorotational instability (MRI), these fluctuating parts should be correctly treated in the general relativistic point 
of view which may require the extended causal thermodynamics such as the Israel-Stewart theory ( Israel fc Stewart 19791 ) 
where the cau sality violating infini te signal speeds are eliminated. The hydrodynamical equations based on such theory are 
formulated by Peitz fc Appl ( 19981 ). If the extended causal thermodynamics for the magnetohydrodynamical fiows can be 
used, the problems of the causality of the viscous fiows at the event horizon will be clearly resolved in the future. 

The another limitation of the viscosity prescription used in this paper is shown in Fig. 1101 In this figure, we show the 
transonic solutions of ADAF calculated by using the Kerr-Schild coordinate {silid lines) and the Boyer-Lindquist coordinate 
[dashed lines) for a/m = 0.95 and 0.99999. The parameters are same as Fig. [8l For the case of a/m — 0.95, the two solutions 
have almost same results. For the cases of lower spin parameters than a/m — 0.95, the same solutions by using the two 
coordinate are obtained. However, in the case of a/m — 0.99999 shown in Fig. 1101 the two results are not exactly same. This 
feature show that the viscosity prescription used in this paper is not perfectly coordinate invariant. When we introduce the 
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causality limited viscosity, we use some special frame in order to evaluate the physical quantities such as the shear stress S or 
the factor fc. However, there is no guarantee that these physical quantities introduced in the specific reference frame have the 
invariant feature with respect to the coordinate transformation or the frame transformation. If we can define and introduce 
the invariant viscosity, the accretion flows calculated by the two procedures produce the same results such as the ideal flows. 
So, the viscous flows presented in the previous sections are not very exactly reproduced from the calculations based on the 
Boyer-Lindquist coordinate such as the past studies. The invariant viscosity prescription will be made by the extended causal 
thermodynamics stated above. 

From here, we give the conclusions of the present study. In the present study, we give the basic equations and the 
calculation method for the horizon-penetrating transonic accretion disks or flows in the equatorial plane written by the Kerr- 
Schild coordinate where there is no coordinate singularity at the event horizon. Based on these formalism, we calculate the 
transonic solutions of these types of the accretion flow models from the outer region to inside region of the event horizon; the 
ideal isothermal flows, the ideal and the viscous polytropic flows, and the advection dominated accretion flows (ADAFs) with 
the relativistic equation of state, the adiabatic accretion disks, the standard accretion disks and the supercritical accretion 
disks. In this study, we use two types of the causal viscosity prescriptions. One is the simple treatment of the kinematic 
viscosity and the other is based on the shear stress measure in the FRF. When we use the causal viscosity prescription 
based on the shear stress measured in the FRF, the boundary condition for the transonic accretion flows is also given at the 
viscous point where the accreting radial velocity is nearly equal t o the viscous diffusion velocity. By using the causal viscosity 
prescription, the ADAF transonic solution s are firstly obtai ned by Gammie fc Pophaml ( 19981 ') for the general relativistic flows 
and later for the pseudo-Newtonian flows ( Takahashi 2007 ) . 

Based on the solutions obtained in the present study, we calculate the physical values for the transonic solutions of the 
these disks around the rotating black hole just on the event horizon and inside the horizon. These solutions are obtained 
for both non-rotating and rotating black holes. In general, the accretion flows calculated by using the Kerr-Schild coordinate 
plunge into black hole with finite three velocity smaller than the speed of light even at the event horizon or inside the horizon, 
and the angular velocities at the horizon are different from the angular velocity of the frame-dragging due to the black hole's 
rotation. These features are different from the results obtained by using the Boyer-Lindquist coordinate with the coordinate 
singularity at the horizon. 

By using the formalism presented in the present study and adding the required physics, we can basically calculate the 
another types of more realistic accretion flows including the radiatively inefficient accretion flows (RIAF) in the galactic 
center, the supercritical accretion disks which is sometimes assumed in the center of the black-hole X-ray binaries or Syfert 
galaxies, and the hypercritical accretion flows or the neutrino-dominated accretion flows (NDAFs) in gamma-ray burst. Also, 
by using the formulations and the calculation methods in this study, the accretion flows inside the black hole satisfying the 
boundary conditions outside the black hole can be calculated. Although the accretion flow structure inside the event horizon 
can not be seen directly by the observer outside the event horizon, by combining the constraints obtained by the future 
observations for the regions just outside the event horizon of the black hole candidate such as the massive black hole at the 
galactic center with the theoretical calculations such as this study or more sophysticated study inside the black hole, we can 
know the accretion flow structure inside the event horizon of the black hole in the real world in the future. Especially, by the 
near-future observations by the radio interferometer such as e.g. VSOP-2 for radio (Hirabayashi et al. 2005), MAXIM for soft 
X-ray (see MAXIM web page: http://maxim.gsfc.nasa.gov/ ), the direct-mapping of the black hole shadow in the RIAF in the 
galactic center will be performed (e.g. Falcke et al. 2000; Melia fc Falcke 2001 for review, Melia 2003a, 2003b and references 
therein) and give the information of the strong-gravity region as the resolved images around the shadow. It is known that 
such images will give the physical information of the black hole itself or the accretion flows in the strong-gravity region (e.g. 
Cunningham & Bardeen 1972, 1973, Bardeen 1973, Takahashi 2004, 2005, Broderick fc Loeb 2005, 2006, Broderick fc Narayan 
2006, Zakharov et al. 2005, Yuan et al. 2006) . However, so far, the image of the black hole shadows calculated by using the 
general relativistic transonic flows of the RIAF have not been performed. Our calculations presented in this study can be also 
applied to such calculations with the radiation mechanisms and the required physics. 
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APPENDIX A: METRIC COMPONENTS 

The components of gfj,„ and are calculated as 



and 



9-=( ^'ii ./^iS.. ). (A2) 



The explicit forms of the non-zero components of g^i/ and gr^" are calculated as 

/ 2mr\ 2mr 2marsm^0 2mr 

gtt = - I 1 2~ j ' 9tr=grt=—^, 9t4>=9<t>t = ^ 1 SrT- = l + -^, 

. 2 2mr \ ^ A sin^ 9 , . „, 

gr4' = g<t,r = -asm e [1 + \ , gee = T,, g^^ = — - — , (A3) 

and 

where nonzero components of /J, and 7'^ arc 

2mr „ 2marsin^6l A rd, sr a ee I a,a, 1 

= 5*. = ^, /3,=5., = g , 7 = s(s + 2mr) ' = = E' ^ = S' ^ = 

Here, we use Pi = fijlS^ and 7**7fcj = <5] where 5j represents the Kronekker delta. 



APPENDIX B: CONGRUENCES FOR THE OBSERVER DRAGGING WITH THE BLACK HOLE'S 
ROTATION 

Here, we show that the congruences for the observer rotating with the angular velocity of the frame-dragging written by the 
Kerr-Schild coordinate have the singularity at the event horizon. We consider the observer rotating with the angular velocity, 
u! = 2mar/A, which is the angular velocity of the frame dragging due to the black hole's rotation. For such observer, by using 
the normalization condition u'^u^ = —1, the contravariant components of the four velocity are described as 

u*=(^^Y\ u''=u'=0, u^=u:u\ (Bl) 

Since A < 0, E > and A > ^ within the event horizon, the component w* become imaginary. Then, in this study, the 
congruences for the observer moving with are not used when the transformation of the physical quantities between the KSF 
and the FRF. 



APPENDIX C: TRANSFORMATION BETWEEN THE KERR-SCHILD FRAME AND THE FLUID'S 
REST FRAME BY TETRADS 



First, we give the tetrad components connecting the KSF and the LNRF made by the congruences of the observer with 
= — a(5^. For such congruences, we have u* = oT^ and u*' = —aT^ji^ (k = r, 6, <f). According with such congruences, the 
metric can be expressed as 

2 

(CI) 



ds^ = -a^dr + — (dr + P^'dty + -leedO" + 7^^ 



,+ Ill(d,. + /3'-dt) 

700 



( ^ ) 



dr + 



E(E + 2mr) 



dr + 



2mr 
T. + 2mr 



dt + Ed6»^ + 



^sin^6i 



ujdt- ^(E + 2mr)dr 



(C2) 
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where we use 7'"'' = 7</>0/(7rr7,#></> — 7r,;,)- By using the tetrad transformation ^ where the hat denote the components 

measured in the LNRF having the orthonormal tetrad basis, the metric can be calculated as 



ds = — 



+ 



'dx" 



(C3) 



Then, the components of the tetrad e„ " connecting between the KSF and the LNRF are calculated as, 
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(C4) 
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By using e'"' = g'^^e^ " and - = rj.^e'^^ where (?;/ii>) is the metric of the Lorentz frame, i.e. diag(?7^j>) = (—1, 1, 1, 1) and 
the non-diagonal components of r]jj,c, are null, we also have the components for the tetrad e** . calculated as 
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Since we now consider the accretion flows in the equatorial plane, we assume that the physical quantities measured 

in the LNRF arc transformed to the physical quantities measured in the FRF by two-dimensional Lorentz transformation 
with the radial velocity Vr and the azimuthal velocity t),/,. Here, we assume that the FRF moves with the radial velocity Vr 
and the azimuthal velocity with respect to the LNRF. Inversely, the LNRF moves with the radial velocity —Vr and the 



azimuthal velocity —v^ with respect to the FRF. The tetrads ' and e"(x^ connecting between the LNRF and the FRF are 
two-dimensional Lorentz transformation with the radial velocity Vr and the azimuthal velocity v,^ measured in the LNRF. 
The transformation matrix c"{X) are described as 
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where 7 = (1 
calculated as 
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'^^"^ . Since both the LNRF and the FRF are orthonormal, the transformation matrix e- are 
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Since now we have both the tetrads between the KSF and the LNRF and the tetrad between the LNRF and the KSF, the 



tetrad connecting the KSF and the FRF, e.g. and e^^ 



, are calculated as 
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e'*(,)=e'',e\.), ^ = e/e. . (CIO) 
By using these tetrad, we can transform the physical quantities between the KSF and the FRF as 

and so on. Lowering and raising the index /i of the tetrads are done by g^n, and g'"^. On the other hand, lowering and raising 
the index /x of the tetrads are done by and j;*'^''''' which are the metric of the Lorentz frame. 



APPENDIX D: VELOCITY FIELDS 



By using the tetrads derived in Appendix[C] in this appendix, we first derive the velocity fields u'^ and in the KSF described 
by the radial velocity Vr and the azimuthal velocity i)0 of the FRF with respect to the LNRF. Since m^*' = — M(t) = 1 and 
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= U(fc) = (fc = r. 
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the four velocity is calculated as u'^ = e** j.^jit'"' = e'^ ^^^u^ ' = 
Then, we have 
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From u*. It'' and '"^(^ we can derive Eq. ([7]) as 
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where we use 6 = -k /2. 
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APPENDIX E: SHEAR RATE 

In this appendix, we calculate the coefficients Ur, Ou and at which gives the shear rate cr(= (y(^r){<f>)) measured in the FRF by 

du'' M 

a — ar + (Ju— \- (Te — . 

ar ar 

The shear rate a(r){4,) in the FRF is calculated from the shear rate a^j^^ in the KSF as cr{r)(c/)) = e^(^)e'^(<^)Crfji^- Here, non-zero 
components of the tetrads are e^'^^.j = t, r, (j)) and e'''(^j (v = t, r, 0) which can be calculated by the tetrad connecting 

the KSF and the LNRF, e*"- , and the tetrad connecting the LNRF and the FRF, e\a) ~ ^) 4') given in Appendix [Cl Since 
o'fj.v = cTi/ji, we need to calculate six components of the shear rate, i.e. att, (^tr, otcp, (Jrr, (^r<t, and a^^. The shear rate ct^^ is 
calculated as the traceless part of the deformation tensor which is calculated as cr^^ = {u^-ah'j^ + Uv-ah'^)/'i — Qh^v/?>- We 
give the covariant derivative for and the four acceleration a^j which are directly used for the calculations of the shear rate 
(jfj^v The non-zero components of are calculated as 

Uf.t ^ \gtt.rU , Ur-t ^ —\{gtt.r + ^gt4,.r)u , U^-t ^ \gt4,,rU' , Ut;r = — ^ ~ {gtt,r + ^gtci>,r)u\ 

2 2 2 ar 2 

1 r / , \ t dl 1 , \ t 1 r 

Ur;r — -^Qrr.rU — (gtr.r + i.lgr4,,r)U , lt0;r = ^ ~ 2^9t<t>,r + ilg^4,,r)U , UQ-B — -^ge0,rU , 

'Ut;<i, — -gt<j>,rU' , Ur;<jy — ~-{gt<j>,r + ^g4,<t>,r)u , U^,^, = -g4,4,,ru' . (El) 

The covariant components of the four acceleration, = u^^.^u'^ , are calculated as 

at = — m''^, ar ^ u''^^ — ^grr,r{u'')'^ ~ {gtr,r + gr4,,r^)u''u' — ^(u)'^g^^^ri^ — ^i){^ — ^K), «fl = 0, a0=u'"^.(E2) 
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We next calculate the expansion O = ti^j-y. The covariant derivatives and u''.^ are calculated as 

1 1 

w*;t = -^(9"gtt,r + g^'''gt<i,,r)u - -^g^''{gtt,r + gt^.r^)u , (E3) 

;r = 2^ grr,rU + (g gtr,r + 9 ^ gr^,r)U + g (ff",'' + 5t0,rS 2) + 9 (gt,#,,r + ff0^,ri 2) J 1i , (E4) 

;9 = 2^ flee,'-'" , (E5) 



"■^i^ = ■^{g*'^gt4>,r + g''"^g4«t>,r)u'' - -g'''^{gt4,,r + g4,g>,r^)u. (E6) 



■^yg gt4>,r + g g<t><i>,r)u -- 

Then, the expansion is calculated by the form 



e = ^^ + e., (E7) 

dr 

where Or is defined as 

r 

r = — {g gtt.r + g grr.r + g gge,r+g g4,4,,r + 2g gtr.r + 2g gt^,r + 2g gr4,,r) (E8) 

We finally calculate the shear rate o-(r)(<t>) which can be calculated as cr(r)(0) = ""r + <Ju{du^ /dr) + ai{d£/dr). By using the 
tetrads are given in Appendix[Cl the shear rate measured in the FRF is calculated by the shear rate (j,^^ measured in 

KSF as 

cr(^)(0) = e''(^je''(<^)cr^^ = + u^;^)e''(^je''(^) - ^ej''^ e^^^y (E9) 

By substituting ti^;„„ given by Eq. HE1|) into Eq. (|E9|) . we can analytically calculate the coefficients (Tr, (Ju and ae. 



APPENDIX F: CALCULATIONS FOR £^ 
Fl Derivation of £* in Kerr-Schild coordinate 

The explicit form for it, is obtained by the direct calculation of £^ = r'^ {Fg^UiiU'^)! where {Fg,^'u-iiu'^)i is calculated from 
rgj^ii^u" until the order of cos^. 2rg^Ui_iu'^ is calculated as 

'^^ei^'^i^^'' = gtt,B{u)^ + grr,0{u)^ + + +2gtr.euu + 2gt,i,.0uu'^ + 2gr,ii.9U u*' , (Fl) 

where the derivatives of the metric components with respect to 9 are calculated until the order of cos 6 around the equatorial 
plane 6 = n /2 as 



gtt,e = gtr,e = 



2ma 



(2 cos 61), 



gt<t>,e 



2ma 
r 



gr4 



2m 
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and ii', u' and u'* are expressed hy £, £ and 

t c fl , 2m\ 2m 
u = £ 1 H H 



r / 



r 2m A a 

U = 1 + —-Mr + 



= — (^ + attr). 



(F2) 
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By substituting ii', vT and u'* given by Eq. (|F3p . we obtain the equation including the term with [urY which is calculated 
from u^Uu ~ —1 as 



A 



{uy = £(1 + 



2m 



r / 



2Ur 



2m ^ a 



\ r 



From Eq. JfTIi with Eqs. ^2^, (|F3} and ((FH), we can finally obtain (FJ^^m^u")! as 



-a\£^ 



(F4) 



(F5) 



where the coeffic ients of Ur become null. 

As shown in Abramowicz. Lanza fc Percivall ( 1997 ). in the Boyer-Lindquist coordinate we obtain the same expression as 
Eq. (|F5p . This is because the specific energy £ and the angular momentum t have the same expression for both coordinate, i.e. 
iSbl = £ks and ^bl = ^ks, as shown by the transformations of four velocities given by Eq. pip . Here, "BL" and "KS" denote 
the physical quantities calculated by using the Boyer-Lindquist coordinate and the Kerr-Schild coordinate, respectively. 
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F2 Derivative of i?, with respect to r 

We express the derivative dlt/dr by the combination of dvl' /dr and dijdr as 

^=,: + C^+4?, (F6) 
ar ar dr 

where it and ii are calculated as 

ll = -'^Sr, r,=-^£u, ei = -^£e + ^. (F7) 

<-* -c* 

Here, £r, £u and are defined to have the relation d£ / dr = £r + £„{dvr / dr) + £(,{(./ dr) and are calculated fully analytically 
or numerically. 
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